This is a tutorial that shows how to make a map of SARS-CoV-2 RNA sewage sampling results in the Netherlands. We’ll start by importing the data. The data comes from this website (Description only provided in Dutch😢) : https://data.rivm.nl/meta/srv/chi/catalog.search#/metadata/a2960b68-9d3f-4dc3-9485-600570cd52b9?tab=relations
For a more general explanation of SARS-CoV2 sewage sampling in the Netherlands (In English!😊): https://www.rivm.nl/en/covid-19/sewage
#make sure you've got the right working directory!
# getwd() and setwd()
nld_sewageRNA <- read.csv("nld.covid.sewage.csv")
#these are spatial packages we'll need to load
library(sf)
library(GISTools)
library(rgdal)
library(sp)
library(tmap)
#this is how you load the Netherlands shape file
# It's important that all the other .dbr, .prj files are stored in the same directory as the .shp file
nld.sf <- st_read("NLD_adm2.shp")
## Reading layer `NLD_adm2' from data source
## `C:\Users\rshea\Desktop\old computer\Thesis code\NLD_adm2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 491 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3.360782 ymin: 50.75517 xmax: 7.29271 ymax: 53.55458
## Geodetic CRS: WGS 84
#here's what the shape file looks like
plot(nld.sf$geometry)
Next we will convert the coordinates from the excel sheet into a spatial object in the sp package and use the rest of the data from the spreadsheet to make a SpatialPointsData Frame
RNA.coords <- data.frame(x=nld_sewageRNA$X_coordinate, y=nld_sewageRNA$Y_coordinate)
nld_sewageRNA.sp <- SpatialPointsDataFrame(RNA.coords, nld_sewageRNA)
Now we need to assign a projection to the spatial points in order to map. The website that had the sewage sampling data provided an epsg code EPSG: 28992
Which you can then use to look up a corresponding proj4 code for on this website: https://spatialreference.org/ref/epsg/?search=28992&srtext=Search
There’s also this website that provides more info about the process:
Then you can use the proj4string(object)<- “+[your proj4string]” command to set the correct projection
proj4string(nld_sewageRNA.sp)<- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs "
## Warning in showSRID(uprojargs, format = "PROJ",
## multiline = "NO", prefer_proj = prefer_proj): Discarded
## datum Unknown based on Bessel 1841 ellipsoid in Proj4
## definition
#then we can convert the sp nld_sewageRNA.sp points to simple features (sf)
#I'm not sure if this step is actually necessary...,but generally speaking sf is easier to work with....but unfortunately also newer so some earlier spatial methods haven't yet been adapted to work in sf, so you still have to run things in sp
nld_sewageRNA.sf <- st_as_sf(nld_sewageRNA.sp)
For the sewage sampling stations mapping we use thematic mapping and we set the mapping mode to view instead of plot. That’s the step that makes an interactive webmap instead of just a static picture
This new map makes the following changes
Dots are proportional in size to the amount of RNA in the sample
The color of the dot corresponds to the date of the sample
The bubbles were made semi-transparent using the alpha term
#this creates a separate layer for names which gives you more choices for things to turn and off in the viewer
nld.admin.labels <- nld.sf
#if this is changed to "plot" you get a static map
tmap_mode("view")
## tmap mode set to interactive viewing
#this creates the admin regions
tm_shape(nld.sf)+
tm_polygons(col="purple", alpha = 0.05)+
#alpha controls the transparency of the polygons
#this adds the sampling result points
tm_shape(nld_sewageRNA.sf)+
tm_bubbles("RNA_per_ml",col="Date_measurement", alpha=0.5)+
#this creates the names layer
tm_shape(nld.admin.labels)+
tm_text("NAME_2", size=1)
## Warning: Number of levels of the variable
## "Date_measurement" is 170, which is larger than
## max.categories (which is 30), so levels are combined.
## Set tmap_options(max.categories = 170) in the layer
## function to show all levels.
## Legend for symbol sizes not available in view mode.
#if this is changed to "plot" you get a static map
tmap_mode("view")
## tmap mode set to interactive viewing
#this creates the admin regions
sewageRNA.map <- tm_shape(nld.sf)+
tm_polygons(col="purple", alpha = 0.05)+
#alpha controls the transparency of the polygons
#this adds the sewage sampling points
tm_shape(nld_sewageRNA.sf)+
tm_bubbles("RNA_per_ml",
col="Date_measurement",
palette= "PuBuGn",
alpha=0.5)+
#this creates the names layer
tm_shape(nld.admin.labels)+
tm_text("NAME_2", size=1)
sewageRNA.map
## Warning: Number of levels of the variable
## "Date_measurement" is 170, which is larger than
## max.categories (which is 30), so levels are combined.
## Set tmap_options(max.categories = 170) in the layer
## function to show all levels.
## Legend for symbol sizes not available in view mode.
Here are some useful documentation files:
https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tmap_animation
https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tm_facets
This step merges all the boundaries in the shape file together. We’ll need this file for making our GIF
netherlands.outline.sf <- st_union(nld.sf)
# the code is set to not run to save space in the tutorial, change eval to TRUE or copy and paste this code into R to make the gif
tmap_mode("plot")
netherlands.outline.sf <- st_union(nld.sf)
Nld_sewageRNA.gif <- tm_shape(netherlands.outline.sf)+
tm_polygons(col="purple", alpha = 0.05)+
tm_layout(title = "SARS-CoV2 RNA in sewage at various sites in the Netherlands")+
#alpha controls the transparency of the polygons
#this adds the sewage sampling points
tm_shape(nld_sewageRNA.sf)+
tm_bubbles("RNA_per_ml",
col="red",
alpha=0.8)+
tm_facets(along = "Date_measurement" , sync = FALSE, as.layers = TRUE)
library(gifski)
tmap_animation(Nld_sewageRNA.gif, filename="Nld_sewageRNA.gif", delay = 100, restart.delay = 300 )
Finally save the file as “NLD COVID19 sewage sampling.Rmd” in your working directory. Then load the markdown library using library(rmarkdown) then use the command render(“NLD COVID19 sewage sampling.Rmd”) to create the HTML page in your library
This code can be used to save a full screen version of the map: